library(Matrix)
library(aroma.light)
## Registered S3 method overwritten by 'R.oo':
##   method        from       
##   throw.default R.methodsS3
## aroma.light v3.14.0 (2018-09-04) successfully loaded. See ?aroma.light for help.
countsFull <- readMM("~/Downloads/GSE128889_RAW/GSM3717977_SCmurinep12_matrix.mtx.gz")
ft <- read.table("~/Downloads/GSE128889_RAW/GSM3717977_SCmurinep12_genes.tsv.gz")
bc <- read.table("~/Downloads/GSE128889_RAW/GSM3717977_SCmurinep12_barcodes.tsv.gz")
rownames(countsFull) <- as.character(ft[,2])

Preprocessing, filtering and trajectory inference

keep <- rowSums(countsFull > 1) >= 400
counts <- countsFull[keep,]
rownames(counts) <- ft[keep,2]
countsFQ <- normalizeQuantileRank(as.matrix(counts))
rownames(countsFQ) <- ft[keep,2]

pca <- prcomp(t(log1p(countsFQ)), scale. = FALSE)
rd <- pca$x[,1:8]

plot(rd[,1:2], pch=16, cex=1/2)

plotByGene <- function(rd, geneCount, ng=10, main=NULL, ...){
  pal <- wesanderson::wes_palette("Zissou1", n=ng, type="continuous")
  gg <- Hmisc::cut2(geneCount, g=ng)
  plot(rd, pch=16, cex=1/2, col=pal[gg], main=main, ...)
}

# progenitors
plotByGene(rd[, 1:2], countsFQ["Dpp4",], ng=2, main="Dpp4")

plotByGene(rd[, 1:2], countsFQ["Wnt2",], ng=4, main="Wnt2")

# group 2 cells
plotByGene(rd[, 1:2], countsFQ["Icam1",], ng=4, main="Icam1")

plotByGene(rd[, 1:2], countsFQ["Dlk1",], ng=4, main="Pref1") #Pref1 is also called Dlk1

# group 3 cells
plotByGene(rd[, 1:2], countsFQ["Clec11a",], ng=4, main="Clec11a")

# group 4 cells
plotByGene(rd[, 1:2], countsFull["Wnt6",], ng=2, main="Wnt6")

# cluster
set.seed(9)
cl <- kmeans(rd, centers = 10)
color = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
pal <- sample(color, 12)
plot(rd[,1:2], pch=16, cex=1/2, col=pal[cl$cluster])
legend("bottomleft", legend=1:10, pch=16, col=pal, bty='n')

# cluster 3 and 7 are progenitor cells
# cluster 1, 5 and 10 are group 2 cells
# cluster 2 and 6 are group 3 cells
keepCells <- which(cl$cluster %in% c(1,2,3,5,7,10))
countsFQ2 <- countsFQ[,keepCells]
rd2 <- prcomp(t(log1p(countsFQ2)), scale. = FALSE)

### umap
library(umap)
rdUmap <- umap(rd2$x[,1:20])
plot(rdUmap$layout, pch=16, cex=1/3)

rafalib::mypar(mfrow=c(3,2))
# progenitors
plotByGene(rdUmap$layout, countsFQ2["Dpp4",], ng=2, main="Prog: Dpp4")
plotByGene(rdUmap$layout, countsFQ2["Wnt2",], ng=4, main="Prog: Wnt2")

# group 2 cells
plotByGene(rdUmap$layout, countsFQ2["Icam1",], ng=4, main="G2: Icam1")
plotByGene(rdUmap$layout, countsFQ2["Dlk1",], ng=4, main="G2: Pref1") #Pref1 is also called Dlk1

# group 3 cells
plotByGene(rdUmap$layout, countsFQ2["Clec11a",], ng=4, main="G3: Clec11a")

# group 4 cells
plotByGene(rdUmap$layout, countsFull["Wnt6",keepCells], ng=2, main="G4: Wnt6")

# cluster
set.seed(12)
pal <- wesanderson::wes_palette("Darjeeling1", n=6, type="continuous")
cl <- kmeans(rdUmap$layout, centers = 6)
plot(rdUmap$layout, pch=16, cex=1/2, col=pal[cl$cluster])
legend("topleft", legend=1:6, pch=16, col=pal[1:6], bty='n')

# slingshot
library(slingshot)
## Loading required package: princurve
lin <- getLineages(rdUmap$layout, clusterLabels=cl$cluster, 
                   start.clus=1, end.clus=c(5,2))
## Using full covariance matrix
plot(rdUmap$layout, pch=16, cex=1/2, col=pal[cl$cluster])
lines(lin, lwd=2)
crv <- getCurves(lin)
## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values

## Warning in regularize.values(x, y, ties, missing(ties)): collapsing to
## unique 'x' values
plot(rdUmap$layout, pch=16, cex=1/2, col=pal[cl$cluster],
     xlab="UMAP1", ylab="UMAP2", bty='l')
lines(crv, lwd=2, col="black")
saveRDS(crv, file="crvUmap.rds")

Evaluate optimal number of knots

8 knots seem appropriate.

library(tradeSeq)
## tradeSeq has been updated to accommodate singleCellExperiment objects as output, making it much more memory efficient. Please check the news file and the updated vignette for details.
set.seed(3)
aicMat1 <- evaluateK(countsFQ2, k=3:10, nGenes=200, sds=crv)

set.seed(4)
aicMat2 <- evaluateK(countsFQ2, k=3:10, nGenes=200, sds=crv)

Model fitting and inference

sce <- fitGAM(countsFQ2, sds=crv, nknots=8)
saveRDS(sce, file="~/data/sceAdipose.rda")

Progenitor cell markers other than Dpp4 and Wnt2 (startVsEndTest)

seRes <- startVsEndTest(sce)
ose <- order(seRes$waldStat, decreasing=TRUE)
which(rownames(seRes)[ose] == "Wnt2") # ranked 263
## [1] 232
which(rownames(seRes)[ose] == "Dpp4") # ranked 450
## [1] 449
head(seRes[ose,], n=16)
##          waldStat df pvalue
## Igfbp7  2360.6550  2      0
## Mfap5   1718.0559  2      0
## Clec3b  1563.6112  2      0
## Pi16    1557.9152  2      0
## Apoe    1324.7050  2      0
## Akr1c18 1296.7321  2      0
## Igfbp6  1196.7044  2      0
## Fn1     1181.2875  2      0
## Anxa1   1118.8296  2      0
## Fbn1    1107.9533  2      0
## Bgn     1098.3094  2      0
## Ly6a    1023.9972  2      0
## Gap43    985.6586  2      0
## Fstl1    941.2143  2      0
## Ly6c1    875.0512  2      0
## Pcolce2  777.1175  2      0
rafalib::mypar(mfrow=c(3,2))
sapply(1:6, function(ii){
  gene <- rownames(seRes)[ose][ii]
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL

## NULL
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
sapply(7:12, function(ii){
  gene <- rownames(seRes)[ose][ii]
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL

## NULL
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
#figure for paper
rafalib::mypar(mfrow=c(3,2))
genes <- c("Dpp4", "Wnt2", 
           "Pi16", "Akr1c18",
           "Fn1", "Fbn1")
sapply(genes, function(gene){
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene, bty='n',
                   xlab="UMAP1", ylab="UMAP2"))
})
## NULL
## NULL
## NULL
## NULL
## NULL

## NULL
## $Dpp4
## NULL
## 
## $Wnt2
## NULL
## 
## $Pi16
## NULL
## 
## $Akr1c18
## NULL
## 
## $Fn1
## NULL
## 
## $Fbn1
## NULL

Differentiated group 2 and group 3 markers (diffEndTest)

Note that this is only restricted to the end points and is therefore limited.

deRes <- diffEndTest(sce)
ode <- order(deRes$waldStat, decreasing=TRUE)


head(seRes[ode,], n=10)
##          waldStat df pvalue
## Meox2    388.2344  2      0
## Igfbp7  2360.6550  2      0
## Cst3     525.6228  2      0
## Nr4a1    203.5429  2      0
## Igfbp6  1196.7044  2      0
## Tmsb4x   712.0800  2      0
## Col14a1  702.0305  2      0
## Cpe      247.7400  2      0
## Cd81     244.5232  2      0
## Bgn     1098.3094  2      0
rafalib::mypar(mfrow=c(3,2))
sapply(1:6, function(ii){
  gene <- rownames(deRes)[ode][ii]
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL

## NULL
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL

earlyDETest for entire cluster of group 2 and 3 cells

plotGeneCount(crv, countsFQ2, models=sce, clusters = cl$cluster)

resedt <- earlyDETest(sce, knots=c(3,7))
oedt <- order(resedt$waldStat, decreasing=TRUE)

head(resedt[oedt,], n=10)
##         waldStat df pvalue
## Dlk1    3521.369  7      0
## H19     3028.316  7      0
## Itm2a   2775.809  7      0
## Mgp     2524.158  7      0
## Rarres2 2436.242  7      0
## Meox2   2188.965  7      0
## Postn   2085.855  7      0
## Fabp4   2062.955  7      0
## Cst3    1934.916  7      0
## Igfbp6  1766.272  7      0
rafalib::mypar(mfrow=c(3,2))
sapply(1:6, function(ii){
  gene <- rownames(resedt)[oedt][ii]
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene))
})
## NULL
## NULL
## NULL
## NULL
## NULL

## NULL
## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
#figure for paper
rafalib::mypar(mfrow=c(1,2))
genes <- c("Mgp", "Meox2")
sapply(genes, function(gene){
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene, bty='n',
                   xlab="UMAP1", ylab="UMAP2"))
})
## NULL

## NULL
## $Mgp
## NULL
## 
## $Meox2
## NULL
#figure for paper
rafalib::mypar(mfrow=c(1,2))
genes <- c("H19", "Col14a1")
sapply(genes, function(gene){
  print(plotByGene(rdUmap$layout, countsFQ2[gene,], ng=4, main=gene, bty='n',
                   xlab="UMAP1", ylab="UMAP2"))
})
## NULL

## NULL
## $H19
## NULL
## 
## $Col14a1
## NULL